Low-lying bifurcations in cavity quantum electrodynamics 



Michael A. ArmerQ and Hideo Mabuchi 
Physical Measurement and Control 266-33, California Institute of Technology, Pasadena, California 91125 USA 

(Dated: February 1, 2008) 

The interplay of quantum fluctuations with nonlinear dynamics is a central topic in the study of 
open quantum systems, connected to fundamental issues (such as decoherence and the quantum- 
classical transition) and practical applications (such as coherent information processing and the 
development of mesoscopic sensors/amplifiers). With this context in mind, we here present a com- 
putational study of some elementary bifurcations that occur in a driven and damped cavity quantum 
electrodynamics (cavity QED) model at low intracavity photon number. In particular, we utilize the 
single-atom cavity QED Master Equation and associated Stochastic Schrodinger Equations to char- 
acterize the equilibrium distribution and dynamical behavior of the quantized intracavity optical 
field in parameter regimes near points in the semiclassical (mean-field, Maxwell-Bloch) bifurcation 
set. Our numerical results show that the semiclassical limit sets are qualitatively preserved in 
the quantum stationary states, although quantum fluctuations apparently induce phase diffusion 
within periodic orbits and stochastic transitions between attractors. We restrict our attention to an 
experimentally realistic parameter regime. 

PACS numbers: 42.65.Pc,42.50.Lc,42.50.Pq 



I. INTRODUCTION 



Bifurcation analysis is a fundamental aspect of dynamical 
systems theory 0, |2j • It provides a powerful set of tools 
and concepts for the study of bistability, hysteresis and 
related phenomena in natural and engineered systems. In 
some practical applications the theory can be used to en- 
sure operation in a structurally stable parameter range, 
while in others it is used to identify operating points that 
are highly sensitive to specific perturbations. The lat- 
ter objective can arise in scenarios where, for example, 
one seeks to exploit the intrinsic nonlinear dynamics of 
a sensing device to provide amplification. This remark- 
ably robust strategy applies even in nanomechanical 3] 
and single trapped-ion [4j systems, as well as in super- 
conducting circuit implementations of quantum compu- 
tation |5|. However, the role of quantum fluctuations 
in determining minimum noise figures and back-action 
for bifurcation amplifiers is not yet understood. It has 
also long been appreciated that bistability and hystere- 
sis could provide a basis for designing logic devices for 
switching and computation in nonlinear optics [jg; inter- 
est in this subject has been reinvigorated by advances in 
the fabrication of photonic bandgap structures and other 
integrated optical circuits 0, • Implementations based 
on strong coupling to quantum dots may provide access 
to a technologically important regime of attojoule and 
picosecond switching, but such performance would seem 
to imply an energy separation between logical states on 
the order of tens of photons. Hence the role of quantum 
fluctuations again calls for attention. Bifurcation analy- 
sis and design thus have many roles to play in modern 
engineering and applied science, and there is a real need 
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to incorporate stochastic analysis. 

In recent years, there has been growing interest in 
studying bifurcation-like behavior of physical, chemical 
and biological systems that are fundamentally discrete 
and stochastic in nature but that can be well described in 
some relevant limit by nonlinear differential equations or 
maps. In chemical reactions, for example, ordinary dif- 
ferential rate equations (corresponding to the so-called 
law of mass action) provide an accurate model in the 
limit where the number of molecules per species becomes 
infinite (at fixed concentration) Q . As a cell biologist or 
as an engineer designing reaction networks for molecular 
computation, however, one might well be concerned with 
understanding precisely how bifurcation phenomena pre- 
dicted by the rate equations are reflected in the stochas- 
tic dynamics of a small number of reacting molecules. 
Stochastic extensions of rigorous bifurcation theory are 
now being developed by several authors 0, 0] , but ex- 
perimental and numerical investigations of specific sys- 
tems are providing crucial guidance for early develop- 
ment of the field [H [H . 

Stochastic behavior and dynamic nonlinearity natu- 
rally coexist in the quantized setting of cavity quan- 
tum electrodynamics (cavity QED) with strong coupling 
[TH ITt| . While the latter subject may seem a bit esoteric 
in comparison to chemical reaction networks, it has im- 
portant connections to nonlinear-optical signal process- 
ing and offers the additional interest of incorporating 
quantum interference and atom-field entanglement. As 
in the case of chemical reactions, the microscopic equa- 
tions of cavity QED are known to have a continuous and 
deterministic macroscopic limit (corresponding to many 
weakly-coupled atoms), which is related to elementary 
models of the laser. It was recognized even in the early 
days of cavity QED [ijj that progress in the laboratory 
would eventually provide a means of exploring bifurca- 
tion phenomena such as optical bistability in a single- 
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atom and few-photon regime; such experiments are in- 
deed now feasible. An opportunity thus arises for the 
empirical study of bifurcation phenomena in the discrete 
physical limit, providing an intriguing quantum-optical 
counterpart to the systems mentioned above. 

While previous theoretical 0, 0, 0] and experi- 
mental |20| investigations of single-atom bistability have 
largely focused on steady-state observables of the trans- 
mitted o ptic al field, we will here follow the spirit of 
Refs. |2l|,|23 in studying transient signals and stochastic 
jumps observable in the broadband photocurrent gener- 
ated in individual experimental trials. We look in par- 
ticular at a case of absorptive bistability, a supercritical 
Hopf bifurcation, and a subcritical Hopf bifurcation, all 
of which occur with mean intracavity photon numbers of 
order ten. Our principal aims in this paper are to illus- 
trate a systematic approach (building upon Ref. [23L Izjf 1 
to expanding the known inventory of bifurcation-type 
phenomena in single- atom cavity QED, and to highlight 
some conspicuous predictions of the fully quantum model 
as compared to the semiclassical Maxwell-Bloch Equa- 
tions. In so doing we hope to begin to illumine a more 
comprehensive picture of the quantum-classical transi- 
tion in cavity nonlinear optics [25L |26| , bridging what is 
generally known about line ar-g aussian |27j and chaotic 
|28| open quantum systems .29]. 



II. THE DRIVEN AND DAMPED 
JAYNES-CUMMINGS MODEL 

A. Quantum Dynamical Description 

We consider the driven Jaynes-Cummings Hamiltonian 
[3(i| which models the interaction of a single mode of an 
optical cavity having resonant frequency lo c , with a two- 
level atom, comprised of a ground state \g) and an excited 
state |e) separated by a frequency u> a . For an atom-field 
coupling constant go and a drive field amplitude £ , the 
Hamiltonian written in a frame rotating at the drive fre- 
quency oji is given by [h = 1]: 

Ti. = Acd 1 * a+A a (T + <j-+igo(aJ a-— ao + )+i£(a) —a) , (1) 

where A a = uj a — oji and A c = lu c — lui . In equation 
(Q, a is the field annihilation operator and er_ = |<?)(e| is 
the atomic lowering operator. In addition to the coher- 
ent dynamics governed by there are two dissipative 
channels for the system: the atom may spontaneously 
emit into modes other than the preferred cavity mode, 
at a rate 711 , and photons may pass through the cavity 
output coupling mirror, at a rate 2k. Furthermore, we 
model the case of non-radiative dephasing (at rate 7 nr ) 
between the atomic ground and excited states. In the 
analysis to follow, we will be concentrating solely on the 
situation where 7„ r = 0, i.e. purely radiative damping; 
however, j nr is included here to indicate that we are not 
restricted to this case (in particular, the parameteriza- 
tion employed in section ITTT1 will imply a variable dephas- 



ing.) The unconditional master equation describing this 
driven, damped, and dephased evolution is 
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where a z = [er+,er_] measures the population difference 
between the excited and ground states. 

While go measures the coherent coupling rate between 
the atom and the cavity, the rates 711, j nr , and k charac- 
terize processes which tend to inhibit the build up of co- 
herence. The qualitative nature of the dynamics (J2J may 
be determined by two dimensionless parameters which 
measure the relative strengths of the coherent and inco- 
herent processes: the critical photon number 
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and the critical atom number 
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where 7j_ is the transverse relaxation rate given by 7j_ = 
7||/2 + 7, ir . The critical photon number provides a mea- 
sure of the number of photons needed to saturate the re- 
sponse of a single atom. Therefore, in the regime no < 1, 
a single photon inside the resonator may induce a non- 
linear system response. Similarly, the critical atom num- 
ber roughly quantifies the number of atoms required to 
drastically change the resonant properties of the cavity. 
When Nq < 1, a single atom inserted into the cavity 
will have a dramatic effect on the cavity output. The 
so-called "strong coupling regime" of cavity QED, which 
is usually used to denote the regime where the coherent 
coupling dominates over dissipation, is reached when the 
condition (no, No) < 1 holds. 

The master equation @ may be used to find the time 
evolution for any operator acting on the system Hilbcrt 
space. In particular, it will be useful to know the dynam- 
ical equations for (a), (er_), and (cr z ) in order to make 
concrete comparisons with the semi-classical results that 
follow. Using the fact that (O) = Tr[Op] for a system 
operator O, we obtain 

(a) = -K(l + iQ){a)+g {(T-)+£ 
(*'-) = -7x(1 + *A)<ct_)+5o<«t,) (5) 
(**) = -7||((O + l)-250«a t ^_) + <a+a)) , 

with 7_l = 7||/2 + 7„ r , 6 = (w c — lui)/k, and A = (u> a - 
w/)/7±. 

It should be noted that these formulae may be easily 
generalized to the case of N non-interacting atoms each 
coupled to the same mode of the electromagnetic field, 
with coupling constant go- In this case, the Hamiltonian 
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becomes 

N N 

H = A c aW^ A Q (T^(Ti+^ ig {a^ o-°_-ao- i + )+iS{a [ -a) , 
3=1 3=1 

(6) 

and the new master equation is 

p = —i[H,p] + K(2apa'—(vap — pa/a) 

N 

+ 7||/2^(2oip4 - u^alp - pu\ul) (7) 
i=i 

JV 

+ 7™r/2^(crJpcrJ -p) , 

3=1 

where cr^_ is the lowering operator for the jth atom and 

[<7j_, <r_] = Sjk&l- The equations of motion for the oper- 
ator expectations become 

N 

(a) = - K (l + iQ)(a)+g J2(<r-)+£ 

3=1 

(at) = - 7± (l + l A)(^)+g (aai) (8) 
(4) = -7||((^) + l)-25o«aV) + <4a)). 

If we define <7- = Ylj=i a - an d a z = SjLi °z as the 
collective pseudo-spin operators, we arrive at the follow- 
ing set of dynamical equations 

(a) = -n(l + i®){a) +g a (<7-) +£ 

= - 1± (l + iA)(<7_) +g (aa z ) (9) 
K> = -7||((^>+^)-2. 9 o((a t ^-> + ^+a»- 

Therefore, we may think of equation I© as a description 
of the operator expectation dynamics for either a single- 
atom (with N = 1) or multi-atom system. In cither 
case, the coupled equations © are not closed, as they 
contain expectation values of operator products. There- 
fore, we also need the dynamical equations for the higher 
order moments, of which there are an infinite number. 
For purely optical systems, order parameters can often 
be identified so that a system size expansion can yield 
a finite, closed set of equations which are valid in the 
"low- noise" limit (when the order parameter is large). 
Unfortunately, for coupled atom-field systems, there ex- 
ists no suitable choice of system scaling parameters which 
would justify a system size expansion 30]. Furthermore, 
it is known that the quantum fluctuations produced by 
optical bistability can be non-classical even when N 1 
[l5|. and therefore would not fit into the classical mold 
which is the basis of a system size expansion. Never- 
theless it has been demonstrated |3l| that the Maxwell- 
Bloch equations, which will be derived from © below, 
can be brought with some refinements into close agree- 
ment with experiments on absorptive optical bistability 



in a multi-atom system. Indeed said equations are gen- 
erally accepted as a canonical, though somewhat phe- 
nomenological, model for c avity nonlinear optics outside 
the strong coupling regime [2|| |33] • 

B. Semi-Classical Description 

An ad hoc (and somewhat crude) approach obtaining 
a closed set of equations from (0 is to simply factorize 
the operator products, e.g. (aV_) — * (at)(cr_). While 
there is no formal basis for this procedure in general, 
the intuition behind it is that for a large system with 
many weakly-excited atoms, the atom-field correlations 
will tend to zero, allowi ng f or expectations of operator 
products to be factorized [32,|33|]- But it should be noted 
that this approximation is not justified in the case of 
strong driving and certainly not for a single atom. This 
factorization yields 

(a) = -K(l + i@){a) +g {o--) +£ 

= - 7x (l + iA)(a_)+ g {a) (a z ) (10) 

W.) = -7||((^)+iV)-2 ff o«a t )<(7_) + (a + )(a)) 

which are the well known the Maxwell-Bloch equations, 
used to describe the semi-classical evolution of a classi- 
cal field coupled to an atomic medium. The atom-field 
correlations which were discarded in performing the fac- 
torization above will tend to contribute "noise" on top of 
the mean field evolution described by Eq. (|10f) . To put 
these equation into a more common form, we make the 
following definitions: 

x = {a), P=f<0, D = ^{a z ). (11) 
so that (|10f> becomes 

x = -n(l + iG)x + (Ng /2)p + £ 
P = -7x(l + iA)p + 2g xD (12) 
D = -^0 + 1)- g {x*p + p*x) . 

A computationally more practical form of l|12l) , which will 
prove useful in the bifurcation analysis to follow, may be 
obtained by transforming it into a dimensionless set of 
equations. We first make the following change of vari- 
ables 

x^y/riaX, p — y-J — p, D — > —D, (13) 

V 7-L 

followed be a re-scaling of time 

t - t'/ 1± , (14) 

so that we are left with the dimensionless Maxwell-Bloch 
equations [32| : 

x = -k[(l+i@)x + 2Cp-y] 
p = -(1 + iA)p + xD (15) 
D = -~f[D- l + (x*p+p*x)/2] 
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where the complex variables x and p represent the ampli- 
tude of the intra-cavity field and the normalized atomic 
polarization, respectively, D is the (real) atomic popu- 
lation inversion, and y is the amplitude of the exter- 
nal drive field. The cooperativity parameter, C, mea- 
sures the strength of the collective atom-field interaction, 
while k and 7 are, respectively, the cavity field decay and 
atomic spontaneous emission rates, scaled by the atomic 
transverse relaxation rate, "f±: 
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The two detuning parameters O and A are the same as in 
Eq. JSJ • Although there is no way to express the steady- 
state solutions for the dependent variables x, p, and D in 
terms of the independent variables, we can find a simple 
set of equations relating the stationary solutions of the 
problem: 



1 



2C 



e- 



1 + A 2 - 
2CA 



-1 2 



1 + A 2 



1/2 



Pa 



(1 - iA)x s 
1 + A 2 + \x. 



1 + A 2 



1 + A 2 



(17) 
(18) 
(19) 



In Fig. n we plot a typical input-output curve gener- 
ated using Eq. I|17|) . Note that in the range 8.7 < y < 
11.1 the curve displays absorptive bistability, with the 
lower and upper branches of the 'S'-shaped curve sup- 
porting stable solutions (the dashed portion of the curve 
is unstable). 

It is important to note that the above equations de- 
pend upon go and N only through the cooperativity 
parameter, C. Thus, identical behavior is predicted 
for a range of systems with varying atom number and 
go = \/2kj±C/N. Of course, one expects that the 
quantum fluctuations and atom-field correlations that 
are disregarded in the derivation of the Maxwell-Bloch 
equations should begin to matter as N approaches 1. 
A direct comparison of 'system behavior' according to 
(|15fl versus the master equation (J2J with the quantities 
{7, k, C, y, 0, A} held fixed can thus be construed as a 
case study in quantum-(semi)classical correspondence. 
The question of course is exactly what 'system behav- 
ior' should be compared and how; the strategy in what 
follows will be to focus on photocurrent properties near 
bifurcation points of the semiclassical model. We thus 
next discuss a systematic approach to finding interesting 
points in the bifurcation set of the Maxwell-Bloch equa- 
tions, and then review a standard monte carlo approach 
to simulating photocurrents according to the quantum 
model. After presenting some numerical results, we con- 
clude with a discussion of some interesting features of 




FIG. 1: Semiclassical calculation of the intracavity steady- 
state field magnitude \x 3S \ versus drive field y. The dashed 
portion of the curve is unstable. The parameter values are: 
C = 10, k = 0.1,7 = 2, 6 = 0, and A = 0. Arrows indicate 
the evolution of the steady-state solution when the drive field 
y is swept smoothly through a bifurcation point: the state 
originally on the lower (upper) branch moving through the 
bifurcation point y « 11.1 (y ~ 8.7) is attracted to the upper 
(lower) branch. 



the quantum-semiclassical comparison that suggest di- 
rections for further research. 



III. BIFURCATION SET OF THE MEAN-FIELD 
EQUATIONS 

In this section we delineate the process used to find and 
classify bifurcations in the mean-field dynamics described 
by Eq. I|15|) . In particular, in Sec. IIII Al we characterize 
both saddle-node bifurcations and Hopf bifurcations. We 
further differentiate between super- and subcritical Hopf 
bifurcations in Sec. IIII Bl In the former case, the bi- 
furcation will destabilize a (typically) fixed point with a 
local (small amplitude) limit cycle born about the prior 
steady-state. In the latter case, no local limit cycle is cre- 
ated about the destabilized steady-state solution, and the 
system will move to a new (possibly distant) attractor. 
For this reason, subcritical Hopf bifurcations often lead 
to qualitatively more radical results, including regions of 
multi-stability. 



A. Linearization about steady-state 

In order to determine the parameter values that lead to 
bifurcations, we return to equation (|15f) . and linearize the 
system dynamics about steady-state. We consider small 
fluctuations Sx, dp, and SD about steady-state, and set 
x — x ss + 5x, p — p ss + Sp, et cetera. After eliminating 
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terms that are second order in the small fluctuations 

Sx 

J 



/ 5x \ 

Sx* 

Sp 
Sp* 

\sbj 



ox 
Sp 
Sp* 
\SDj 



(20) 



where the Jacobian J is given by 
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/k(l + iG) 2Ck 
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-D ss 1-iA -x 
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The associated characteristic equation will have the form 

A 5 + aiA 4 + a 2 X 3 + a 3 A 2 + a 4 A + a 5 = (22) 
with coefficients given by the following expressions: 
oi = 2 + 7 + 2fc , 

a 2 = fc 2 (l + e 2 ) + (2 7 +l + A 2 +7|:r ss | 2 ) 

+ 2fc( 7 + 2) + AkCD ss , 
a 3 = 7 (l + A 2 + |a; s . s | 2 ) + 2fc(27+l + A 2 +7|x ss | 2 ) 

+ fc 2 (l + 6 2 )( 7 + 2) + 4fcCZ? ss ( 7 + k + 1) 

- lkC(p* ss x ss +p ss x* ss ) , 
a 4 = 2fc 7 (l + A 2 + |a; ss | 2 ) 

+ fc 2 (l + e 2 )(27 + l + A 2 +7|a: ss | 2 ) 

+ 2kCD ss [2k(l - AG) + 2 7 (fc + 1) + 7|x ss | 2 ] 

+ 7 /cC [i(A + k&)(p* s x ss - Pss x* ss ) 

- (k+ l)(p* ss x ss +p..x*,)\ + Ak 2 C 2 D 2 s , 
a 5 = -fk 2 {4C 2 D ss [D ss - (p* s x ss + p ss x* ss )2] 

+ (l + e 2 )(l + A 2 + | 2 ; ss | 2 ) 
+ 4CD SS (1 - A6)} . 

These coefficients may be further simplified by using the 
form of D ss found in (|17|l and the relations 

i( Pss x* ss -p* ss x ss ) = -2A\x ss \ 2 /(l + A 2 + \x ss \ 2 ) 
(p ss x* ss +p* ss x ss ) = 2|a; ss | 2 /(l + A 2 + |a; s;5 | 2 ) , 

so that the a's found in the characteristic equa- 
tion are written explicitly in terms of six parameters: 
C, k, 7,0, A, and |x sa |. 

Needless to say, it is impossible to solve for the eigen- 
values of this system analytically. However, other meth- 
ods that provide analytic tests of stability do exist. Most 
notably, the Routh-Hurwitz criterion provides a set of 
inequalities based on combinations of the a's that can 
be used to determine stability. Unfortunately, this pro- 
cedure is simply too general, and it is ill suited for the 
purpose of determining the boundaries of instability in 
terms of our controllable parameters. We can, however, 



make use of the Routh-Hurwitz criteria to find the fol- 
lowing necessary conditions for stability |34l l35| 

at, a 2 , d3, o 4) a 5 > . 

Furthermore, at a Hopf bifurcation the system must have 
a pair of pure imaginary eigenvalues, Ai,2 = iiw. De- 
manding that the characteristic equation support these 
solutions establishes the following critical condition for a 
Hopf bifurcation 

/ = (a 1 a 2 ~a 3 )(a 3 a i ~a 2 a 5 )~ (aia 4 -a 5 ) 2 = , (24) 

with / > providing another necessary condition for 
stability. 

For the purpose of delineating the instability bound- 
aries, the six inequalities 



ax, 02, a 3 , a 4 , a 5 , / > 



(25) 



are not equally important. Starting from a stable region 
of the parameter space, there are only two ways for the 
steady-state solution to become unstable: (i) a single real 
eigenvalue passes through the origin and becomes posi- 
tive; (ii) a pair of complex conjugate eigenvalues cross 
the imaginary axis (starting from the LHP). For case (i) 
the coefficient a 5 must change signs first, whereas for case 
(ii) it is / that first changes sign. Therefore, if the goal 
is to determine the conditions for a known stable state to 
become unstable, there is no need to consider the other 
necessary conditions and all focus may be placed on 05 
and /. Furthermore, if we are only interested in Hopf Bi- 
furcations, we can also ignore the 05 > condition, which 
determines the boundary for saddle-node bifurcations 
where the steady-state curve displays a turning point 
(this can be seen by noting that d(y 2 )/d(\x ss \ 2 ) cx a 5 , 
so that 05 < indicates bistability.) 

It should be noted again, that the inequalities lt2~5l) 
are only necessary conditions for stability, they are not 
sufficient. For example, the system could have one real 
negative eigenvalue, and two pairs of complex eigen- 
values each with positive real parts, and still have 
01,0-2,03,04,05,/ > 0. However, the stability condi- 
tion / > can be made sufficient for a given region 
in parameter space if we can show that this region is 
'connected' to a known stable region of the space (two 
regions of parameter space are 'connected' if there ex- 
ists a continuous variation of the parameters that moves 
the system from one region, while retaining the sign of 
01, 02, 03, 04, 05, and / through the entire path.) Thus, 
if we know that a particular region of parameter space 
(with ax, ai, 03, 04, as > 0) is connected to a stable re- 
gion, we know that / can serve as a necessary and suf- 
ficient condition for the steady-state solution to undergo 
a Hopf bifurcation. Practically, all this means is that 
starting from a stable state, the first crossing of a surface 
05 = (/ = 0), will drive the system unstable through 
a saddle-node (Hopf) bifurcation. Furthermore, the sta- 
bility condition / > is quite reliable in practice, even 
when we can't show connectedness to a stable region. 
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B. Super- and subcritical Hopf bifurcations 



with 



In order to determine whether a Hopf bifurcation is 
super- or subcritical, the eigenvalues and eigenvectors 
about the bifurcation point must first be found. Among 
the possible reasons for seeking one or the other kind 
are that supercritical Hopf bifurcations can be used for 
resonant nonlinear amplification of small periodic signals 
[36| , and subcritical Hopf bifurcations are likely to indi- 
cate the presence of limit cycles that coexist with other 
attractors. The latter type of scenario may give rise to 
observable 'quantum jumps' among non-fixed point at- 
tractors, which would be an interesting generalization of 
the predictions of Refs. pll l37| . We thus believe that 
the theory for distinguishing super- and subcritical Hopf 
bifurcations merits an extended discussion. Note that 
our expressions below and in the Appendix correct some 
apparent misprints in Ref. j23| . with minor changes of 
notation. 

At a Hopf bifurcation, the linearized system l|21|l has a 
pair of pure imaginary eigenvalues A1.2 = ±ilu, with the 
frequency uj determined by 

= aiQ4 ' a5 . (26) 
aia 2 — a 3 

Thus, the characteristic equation can be factored as 



(X 2 +uj 2 ){X 3 + ai X 2 + b 2 X + b 3 ) =0 



(27) 



62 = fl4/^ 2 , and 63 = a^/ui 2 . (28) 

Solving for the other three eigenvalues yields 

A3 = vw — p/(3vw) — ai/3 

A 4 = v*w-p/(3v*w)-a 1 /3 (29) 

A5 = w — p/(Zw) — oi/3 



where the variables 



v = 
w = 



(-l + »V3)/2 



(30) 



q = -2af/27 + ai6 2 /3- 63 

f>2 



p = -a 2 /3 + 6 2 



are determined from the solution to the cubic equation 
embedded in l|27|l . Following the approach in [23j (and 
noting several corrections), the system eigenvectors, a- u 
may be found in terms of the Xi by solving the linearized 
dynamics 



Jtti = AjQi 



(31) 



where J is the Jacobian in (|21|l . Expressing the results 
in terms of the A,- , one arrives at 



/ exp (ifc) ^ 

exp (— i(pi) 
- exp (ifa) [fc(l + iO) + Xi] /(2Ck) 
- exp {-ifa) [fc(l - iQ) + Xi] /(2Ck) 
V - exp (ifc) {2CD ss k + (1 + iA + A,) [k(l + iO) + X t ]} /{2Ckx) J 



(32) 



where the phase factor, exp(i^), is chosen to preserve 
symmetry in the components and is given by 

i(j>i = I x{2CD ss k + (l-iA + A.Q[fc(l - jgj + A~]y 
~ y x* {2CD ss k + {1 + iA + Xi)[k{l + iQ) + Xi]}' 

(33) 

The set of eigenvectors, {«i}, define a linear transfor- 
mation of variables such that dynamical equations about 
steady-state contain no linear cross couplings. The old 
variables are related to the new variables through the 
relation 

(5q=5z (34) 

where a = [a± u-z 04 as] is the transformation ma- 
trix comprised of the system eigenvectors, and <5q = 
{Sx, Sx* , Sp, 5p* , SD} T are the fluctuations about 



steady-state. Starting from the dynamical equations for 
the fluctuations about steady-state 

Sx = -k [(1 + iG)6x + 2CSp] 
Sx* = {Sx)* 

Sp = —(I + iA)Sp + SxSD (35) 

+ {SxD ss + x ss SD} 
Sp* = (Sp)* 

5D = -7 [SD + (Sx*Sp + 8p*8x) 12] 

+ {l/2(8x*p ss + x* s Sp + p* ss Sx + 5p*x ss )} , 

the transformation of coordinates z = (3 Sq, with /3 
11 eliminates any linear coupling between variables 

ij = XjZj + fij^SxSD + /3ji5x*SD 

- -i(3 jb (5x*Sp + Sx5p*)/2 . (36) 
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Finally, after utilizing the transformation l|34(l , this equa- 
tion may be expressed in terms of the "diagonalized" co- 
ordinates alone: 

5 

Zj = ^3 z j + X! i^33 a ika 5 i + /3j4a 2 ka 5 i 

k,l=l 

- ~/f3 j5 (a 2k a 3 i + ai k a 4 i)/2] z k z t . (37) 

Having converted the system into the form (|37[) . the 
dynamics about a Hopf bifurcation may be reduced onto 
a center manifold [381 ]: since the system dynamics will 
be dominated by the "slow" variables, Z\ and z 2 , the 
flow of the differential equation may be locally approx- 
imated on the surface generated by Z\ and z 2 , with the 
"fast" variables, 2^=3,4,5, represented by a local graph 
Zj = Aj(z\, z 2 ). Furthermore, the local graph, Aj(z%, z 2 \ 
may be approximated by a power series expansion 

Zj = a 20 (j)zf + a 11 (j)z 1 z 2 + a 02 (j)zj + ... , j = 3,4,5 , 

(38) 

so that the reduced dynamics may be approximated by 

zj sa 2iuja 2Q {j)zl - 2iuja 02 (j)zl + ... , j = 3,4, 5. (39) 

The coefficients in (I38|) are determined by substituting 
equations (|38|) and (|39|l into the exact dynamics l|37|) and 
equating like powers in z"z™. 

With aid of the local graph J3HJl, the dynamics may 
be reduced onto the center coordinates associated with 
eigenvalues having zero real part: 

zi = iujzx+b^i^zf+bui^ziz^bo^^z^+hii^zlz^... 

(40) 

with z 2 = z* near the Hopf bifurcation. Writing z\ = 
u + iv, with {u,v} € R, the reduced dynamics of the 
complex variable Z\ may be expressed by a set of real 
coupled equations 

u = -ujv + F(u,v) , 41 < 
v = luu + G(u, v) ' 

with F and G comprised of terms nonlinear in u and v. 
It can be shown |38||. that there exists a smooth change 
of variables which will put equation l|41|l into the normal 
form 

r = r] 3 r 3 + r) 5 r 5 + ... 

6 = uo + e 2 r 2 + e 4 r 4 + ... ' 

with the stability of the bifurcation governed by the sign 
of 773. For 773 < the bifurcation is supercritical, and 
a stable, small amplitude limit cycle is born about the 
newly destabilized steady-state. For 773 > 0, the bifurca- 
tion is subcritical, and no such small amplitude cycle is 
created. In this case, the system may be thrown far away 
from the steady-state solution, onto cither a limit cycle, a 
different branch of the steady-state curve, or some other 
attractor. 

Remaining details of the calculation of 773 are relegated 
to the Appendix. 



IV. QUANTUM SIGNATURES OF 
BISTABILITY AND LIMIT CYCLES 

In this section we briefly review some computational 
tools that can be used to search for evidence of the semi- 
classical attractors in the quantum model. In the follow- 
ing section we present numerical results in which these 
tools are applied at several interesting points in the semi- 
classical bifurcation set. 



A. Steady-state Q-function 

The Q-function (for the intracavity optical field) [3^ 
can be used to characterize the steady-state behavior of 
the driven atom-cavity system. The benefits of using 
the Q-function, over the other common quasi-probability 
distributions, are twofold. Firstly, the Q-function is 
positive-semidefinite, and is therefore better suited for 
making comparisons to classical probability distributions 
(such as the stationary distribution of a classical model 
with noise). Secondly, the value of the Q-function has a 
very simple interpretation in terms of coherent states \a) 
of the intracavity field, since for any given density oper- 
ator p we have Q(a) = (a\p\a). Thus Q(a) may strictly 
be thought of as a probability, which further lends to the 
utility of the Q-function as a tool for making comparisons 
between the quantum and semiclassical descriptions of 
the equilibrium behavior. 

In practice we can find the steady-state density op- 
erator p ss of the master equation (J2J simply by set- 
ting the terms on its right-hand side to zero and solv- 
ing the resulting algebraic equation. We can then use 
Q S s{&) = ( a \pss\(x)-> where there is an implied partial 
trace over the atomic degrees of freedom. 

One naively expects that the Q-function should be 
bimodal when atom-cavity parameters are chosen in 
a bistable region of the Maxwell-Bloch equations [l7| . 
Likewise, the existence of a limit cycle should give rise 
to a ring-shaped Q-function. Below we will show exam- 
ples of such features, but we will also want a method 
for visualizing any possible coherent dynamics or large 
fluctuations buried within these stationary distributions. 
For example we want to be able to show that a ringlike 
Q-function indicates an intracavity field amplitude that 
oscillates coherently, as opposed to randomly diffusing 
around the circle. In the case of bistability we would 
also like to see the switching timescale and to examine 
the sharpness of the switching events. 

B. Quantum Trajectories 

To visualize dynamics in the quantum model in an 
experimentally-relevant way, we turn to the method of 
quantum trajectory simulations [30I |40| . 

The master equation @ generates predictions of the 
unconditioned state of the atom-cavity system. That 
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is, the solutions p(t) of the master equation represent 
the knowledge we can have of the evolving system state 
without utilizing the information obtainable via real-time 
measurements of the output fields (cavity transmission 
and atomic fluorescence). Hence we may gain further in- 
sight into the dynamics by considering quantum predic- 
tions regarding the photocurrents; fortunately, the same 
theory used to derive the master equation provides a pow- 
erful set of tools for statisticall y-fa ithful sampling of con- 
tinuous measurement records |8?j Eflj , and tells us how 
to interpret them as real-time observations of the intra- 
cavity dynamics |4l| . Here we will use such 'quantum 
trajectory' methods for monte carlo simulations of the 
photocurrent generated by homodyne detection of the 
cavity output field [4(j] . 

For the case of homodyne detection of the output fields, 
the stochastic Schrodinger equation (SSE) governing the 
evolution of the (unnormalized) conditional state vector, 
\ip c ), is given by 

d\ip c ) = —i\H — ina^a — i-f ^(7 + er_) \i/j c ) + 

V2Kae i ^\ip e )dQ 1 + (43) 
v /^yI<7_e*|V c )dQ 2 , 

with the Hamiltonian, TC, given in Eq. Q), and the phase 
factor e 1 ^ 1,2 determining the field quadrature being mea- 
sured (e.g. 4>\^2 = for measurement of the amplitude 
quadrature and 4>x,2 = — tt/2 for the phase quadrature.) 
The measured homodyne photocurrents, //j 0mi 2 = rf ^' 2 , 
which are, respectively, the homodyne photocurrents as- 
sociated with the cavity and atomic decay channels, are 
calculated using 

dQi = V2 ^ )c dt + dW, , (44) 

, In , p ~ i( t>2 _±_ n p i<p2\ 

dQ 2 = /,7/ \ dt + dw * » («) 

where dW\ 2 are independent Wiener increments satisfy- 
ing (dW) = and ((dW) 2 ) — dt. Numerical integration 
of equations (|43l44[l are performed using the stochastic 
integration routine incorporated in the Quantum Optics 
Toolbox for MATLAB. 

In any practical experiment, full measurement of the 
atomic spontaneous emission is not actually feasible as 
this would require a detector covering nearly A-k steradi- 
ans of solid angle. Fortunately, the cavity-output homo- 
dyne photocurrent 4 omi generated by monte carlo inte- 
gration of the above SSE (considered on its own with- 
out any reference to the corresponding Ihom 2 ) is sam- 
pled from the same law as the photocurrent one would 
see in an experiment in which the atomic decay channel 
was not measured at all [43. We therefore make use of 
such photocurrent simulations below in our discussion 
of single-atom bistability and Hopf bifurcations. One 
should appreciate however that the conditional state |-0 C ) 
propagated by the SSE is then merely an internal vari- 
able of the monte carlo simulation, and not something 



that could actually be reconstructed (in a recursive esti- 
mation sense) from just the cavity-output photocurrent. 
The best one could do along the latter lines, without 
assuming high-efficiency observation of the cavity decay 
channel, would be to utilize the corresponding stochastic 
master equation (SME) 40j as an optimal quantum filter 
|4l| . Even in a purely theoretical discussion one would 
like to utilize the SME (as in Ref. (3?|]) to generate not 
only realistic photocurrent samples, but also the condi- 
tional quantum states that one could in principle gener- 
ate from them via recursive filtering. Unfortunately such 
numerical procedures are very computationally intensive. 
As we find that adequate indications of the dynamics un- 
derlying bimodal and ringlike Q-functions are provided 
by the photocurrents alone, we have limited our efforts 
to SSE-based simulations. 



V. NUMERICAL RESULTS 
A. Absorptive bistability 

In Fig. n we plot the steady-state intra-cavity field 
magnitude vs. drive field predicted by the (dimension- 
less) semi-classical equations (ITTIl for the case of purely 
absorptive bistability (0 = A = 0). These parameter 
values correspond to g = 1.41, k — -1,7|| = 2,A a = 
0, A c = 0, in the master equation J5J, and a saturation 
photon number, no = 0.25. In Fig.|21we plot the bimodal 
Q-function obtained from the steady solution to the mas- 
ter equation for a drive field y = 11.3 (£ = 0.565), cho- 
sen such that the integrated probabilities in each mode 
of the Q-function are approximately equal. While this 
Q-function indicates that the quantum dynamics show 
bistable behavior, it is interesting to note that the master 
equation produces this bimodal distribution for a drive 
field where the mean-field equations do not predict bista- 
bility (the lower branch in Fig.^disappears at y w 11.1). 
In fact, the Q-function distributions over most of the 
semi-classically bistable region are not bimodal, and only 
become so for 11 < y < 11.5. 

In Fig. |21 we plot the photocurrent l|44|l correspond- 
ing to a measurement of the amplitude quadrature of the 
cavity output field. As expected, the field localization 
brought about by the continuous homodyne measure- 
ment causes the signal amplitude to switch, at stochastic 
intervals, between values consistent with the two peaks 
in Q(a). 

B. Supercritical Hopf bifurcation 

In Fig. 0| we plot the equilibrium attractors of the 
mean-field dynamics l|15l) for a case where the steady- 
state fixed points predicted in (|17f) undergo supercritical 
Hopf bifurcations. Starting on the lower (upper) branch 
of the steady-state curve, as the drive field y is swept 
through the critical point CP\ (CP2) the fixed point is 



Re(u) 

FIG. 2: Steady-state Q(a) for the parameter values in Fig.Q 
and drive field y = 11.3. 
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FIG. 3: Homodyne photocurrent from quantum trajectory 
simulation for parameter values in Fig. [5] 



destabilized by a small amplitude limit cycle, LC\, which 
grows in amplitude, peaking at y w 2800, until finally re- 
combining with and restabilizing the fixed point at CPi 
(CP\). To represent the oscillatory solution that is born 
out of the bifurcation, we plot the steady-state maximum 
field magnitude for a state localized on the stable limit 
cycle, and denote this as LC\. Thus, the plotted curve 
LC\ essentially represents the amplitude plus mean value 
of the limit cycle. 

The parameter values used in Fig. 0] correspond to 
g = l,K= .01,7|| = 2,A a = 1.25, A c = -6 in Eq. ©, 
and a saturation photon number uq — 0.5. Using these 
values, we compute Q(a) for y — 2800 {£ = 19.8), where 
the limit cycle amplitude is maximal. The result is plot- 
ted in Fig. [3] The ring-like shape of the distribution is 
consistent with oscillation of a coherent state in the in- 
tracavity field. This interpretation is further supported 
by the inset in Fig. where we plot the autocorre- 
lation function G [ y\t) = (Y(t)Y(0)) - (y(0]) 2 , com- 
puted using the quantum regression theorem 39] , where 
Y = ^(a^— a) is the phase quadrature amplitude operator 
of the intracavity field. In addition, Fig. [5] displays the 
coherence time of the steady-state quantum oscillations 



500 1000 1500 2000 2600 3000 3500 4000 4500 5000 5500 

y 

FIG. 4: Semiclassical calculation of the intracavity steady- 
state field magnitude \x ss | versus drive field y. The parameter 
values are: C = 50, k = .01, 7 = 2, 9 = -600, and A = 1.25. 
The upper dashed-dotted curve (LCi Represents the steady- 
state oscillation maximum (steady-state magnitude plus mean 
value) associated with the limit cycle formed when the fixed 
point becomes unstable due to a Hopf bifurcation (at CP\ 
and CP?). The lower dashed curve is unstable. 




FIG. 5: Steady-state Q(a) for the parameter values in Fig.Q] 
and drive field y = 2800. 

over a range of drive fields. The results indicate that the 
coherence times depend strongly on the amplitude of the 
limit cycle, LC\, which is again consistent with the idea 
of an oscillating coherent state for the intracavity field. 

It can be seen clearly from the inset of Fig. O that the 
limit cycle comprises an oscillation of the intracavity field 
at a frequency much higher than k. It should thus be dif- 
ficult to see the oscillation directly in the broadband pho- 
tocurrent generated by amplitude-quadrature homodyne 
detection of the cavity output field. In Fig. however we 
plot several power spectra of photocurrent records gen- 
erated in quantum trajectory simulations. For y = 1000 
(below CP\ in Fig. 0} the spectrum shows little or no 



10 



2.6 



2.4- 




1000 1500 2000 2500 3000 3500 4000 4500 5000 

Drive Field [y] 



FIG. 6: Coherence times calculated from steady-state au- 
tocorrelation function, Gy\t), of the quantum mechanical 
phase quadrature. Parameter values same as in Fig. 2] The 
coherence time is estimated by fitting a damped sinusoid to 
Gy' (t). Inset: simulated Gy (t) for the drive field in Fig.|S] 

sign of a coherent peak, but for y — 2800 we see that ho- 
modyne detection of the field amplitude reveals clear evi- 
dence of the limit-cycle oscillation. This demonstrates at 
least a basic correspondence with the semiclassical pre- 
dictions shown in Fig. At y — 5000 (above CP-i in 
Fig. QJ, however, we see that the quantum model still 
exhibits strong oscillations even though the semiclassical 
model predicts a fixed point solution. This persistence of 
the oscillatory behavior at both higher and lower driving 
fields can also be seen in Fig. 

C. Subcritical Hopf bifurcation 

In Fig.|S]we plot the steady-state solutions for a param- 
eter regime where the mean-field equations predict a sub- 
critical Hopf bifurcation. The solid (dashed) curve cor- 
responds to the stable (unstable) fixed points predicted 
by Eq. ijl"7|) . whereas the attractor LCi (plotted dashed- 
dotted) corresponds to a stable limit cycle. Beginning 
on the upper stable branch of fixed points, as the drive 
field is swept through the critical point CP4, the system 
undergoes a subcritical Hopf bifurcation. In the range 
475 < y < 925, the semiclassical equations predict coex- 
istence of a stable fixed point and limit cycle, which is a 
common signature of subcritical bifurcations. Note that 
at y « 925 the limit cycle LC2 is destabilized but the 
fixed point is not. The two arrowed lines in Fig.[S]do not 
represent solutions to the equations, but simply indicate 
which attractor a destabilized state will seek. 

The parameter values used in Fig. [S] correspond to 
g = 4.47, k = .05,7|| = 2,A a = 2,A C = -2.75, and a 
saturation photon number no = 0.025. The small size 
of no indicates that the qualitative behavior of steady- 
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Frequency [k] 

FIG. 7: Power spectra computed from simulated photocur- 
rents for amplitude-quadrature homodyne detection of the 
cavity output field, using parameter values as in Fig. [I] Low- 
est plot: y — 1000, below the bifurcation point; middle plot: 
y = 2800, where the classical oscillation amplitude is maxi- 
mal; upper plot: y — 5000, where the semiclassical model no 
longer predicts a limit cycle. 
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FIG. 8: Semiclassical calculation of the intracavity steady- 
state field magnitude \x ss | versus drive field y. The parameter 
values are: C = 200, k = .05,7 = 2, 6 = -55, and A = 2. 
The upper dashed-dotted curve (LC2) represents the steady- 
state oscillation maximum (steady-state magnitude plus mean 
value) associated with the limit cycle formed when the fixed 
point becomes unstable due to a Hopf bifurcation. The critical 
point CP3 is supercritical, whereas the point CP4 is subcrit- 
ical. Note the coexistence of a stable fixed point and limit 
cycle in the range 475 < y < 925. The lower dashed curve is 
unstable. 



state solutions to the master equation will be dominated 
by quantum fluctuations (over the dynamics implied by 
the mean field equations.) We see that this interpreta- 
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FIG. 9: Steady-state surface (a) and contour plots (b) of Q(a) 
for the parameter values in Fig. |H1 and drive field y = 1000. 



tion is justified by the plot of Q(a) in Fig. EI at a drive 
field y = 1000, near the amplitude maximum of LC'2- 
The features in the surface plot (Fig. Et 1 ) are hardly pro- 
found. The contour plot in Fig. [SJd is somewhat more 
elucidating, and shows signs of the coexistence of oscilla- 
tory, albeit asymmetric, and fixed coherent states. These 
'blurry' results are not surprising, as the small no, and 
relation (|13l) . imply that, in the quantum case, the struc- 
ture implied by the limit cycle and fixed points in Fig. [S] 
should be heavily affected by fluctuations. 

Fig. suggests, however, that the field localization 
provided by continuous homodyne measurement of the 
phase quadrature can reveal signatures of bistability in 
the photocurrent. In particular, in Fig. 1 10b we plot a 
typical trajectory for a drive field y = 1000, as in Fig. [5] 
While the contrast is marginal, as it would have to be 
given the discussion above, the qualitative appearance of 
this simulated photocurrent record is that of oscillations 
interrupted by brief periods of stationary noise (which 
one could attribute to transient localization on the fixed 
point). Such intermittency can also be seen in Fig. llOt;. 
which corresponds to a drive field, y — 1400, well past 
the semiclassical region of multi-stability. 
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FIG. 10: Homodyne photocurrent from quantum trajectory 
simulations with parameter values in Fig. [8] (a) drive field 
y = 500, near the subcritical bifurcation point; (b) for drive 
field y — 1000, near the oscillation amplitude maximum; (c) 
for drive field y — 1400, beyond the region of semiclassical 
bistability. 



VI. CONCLUSIONS 

In this paper we have considered bifurcation phenom- 
ena as focal points for the investigation of quantum- 
(semi) classical correspondence in cavity nonlinear optics. 
We presented a general approach to the characteriza- 
tion of interesting points in the semiclassical (Maxwell- 
Bloch) bifurcation set, starting from the formal methods 
of Refs. |U HH (but with corrections to some of their 
equations as printed) and incorporating numerical simu- 
lations of the cavity QED master equation and (homo- 
dyne) stochastic Schrodinger equation. This approach 
leads to the prediction of self-oscillation and bifurcation- 
type behavior in an experimentally accessible parameter 
regime for single-atom cavity QED, under driving condi- 
tions in which the mean intracavity photon number is of 
order ten. It is interesting to see that such a wide range 
of input-output characteristics are supported in such a 
small state space. 

The results of our numerical simulations point to a 
number of questions for further study. For example 
one would like to understand, in physical terms, what 
determines the correlation timescales found in Fig. |3| 
(bistability) and Fig. (stable limit cycle). In scenar- 
ios where switching occurs between output-signal char- 
acteristics associated with different semiclassical attrac- 
tors (as in Fig. |3| and Fig. ITTjl) . one would also like to 
know what determines the attractors' relative stability 
and how the switch ing events are initiated (in the sense 
of Refs. 0,H3,|33,E3)- Can we understand the physi- 
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cal dynamics that give rise to limit-cycles in single-atom 
cavity QED (as in Ref. HJ), and why they are desta- 
bilized in subcritical Hopf bifurcations? Such investiga- 
tions could be taken as a starting point for the devel- 
opment of quantum feedback control strategies 0, 23 
to stabilize desired fixed points 0, 23 or limit cycles, 
or for inducing switching between them with minimum 
time or dissipated energy (as would be required for the 
types of optical signal processing applications mentioned 
in the introduction). More generally, one can ask about 
the applicability of ideas from classical bifurcation con- 
trol [43, I43. In this context, single-atom cavity QED 
provides an interesting model system for investigations of 
quantum feedback control far beyond the linear-gaussian 
regime [^IHflj . 

Returning finally to the issue of correspondence, one 
curious detail of our simulation results is that near hys- 
teresis loops, bimodal behavior in the quantum model is 
often seen to persist well above (at higher driving fields 
than) the upper switching points of the Maxwell-Bloch 
equations, where saddle-node or subcritical bifurcations 
occur. We have also seen a case (supercritical Hopf) in 
which oscillatory behavior occurs over a wider parameter 
range than is predicted by the semiclassical model. It is 
certainly possible that these effects may be explainable 
by analogy with noise-induced 'postponement' [5ll l52|. 
'advancement' or 'precursors' [53] in classical nonlinear 
systems, or with some other known effect in classical ran- 
dom dynamics [3 EE EE 113- If this is the case, it will 
be interesting to see what clues such analogies provide to- 
wards the development of physically-motivated stochas- 
tic extensions of the Maxwell-Bloch equations that can 
capture global (in the dynamical systems sense) effects 
of quantum fluctuations (thus going beyond what is pos- 
sible with local- linearization approaches, as in Ref. [58j). 
In any case it will be natural to try to relate the machin- 
ery of 'P-bifurcation' analysis [IE EH to our steady-state 
Q-functions. Of course if sensible analogies with post- 
ponement in classical noise-driven systems cannot be es- 
tablished, it will be tempting to ask whether coherence 
or atom-field entanglement play any significant role in 
these or any other cavity QED bifurcation phenomena. 



In any case, we certainly expect that further study of 
bifurcation phenomena in single-atom cavity QED will 
improve our understanding of the way that quantum fluc- 
tuations 'enrich' the mean-field phase portrait of coherent 
nonlinear dynamical systems. 
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APPENDIX A: CALCULATION OF 773 

The critical parameter rj 3 can be expressed explicitly 
[38l E3 m terms of the coefficients in l|4Tj|) and the bifur- 
cation frequency to 

V3 = Re[&2i(l)] - -Im[Ml)Ml)] • (Al) 

U! 

As the calculations for the relevant coefficients are rather 
tedious and time consuming, they are included here for 
posterity: 

oaoC?) = Mj)/(2^ - A,) , 1 = 3 4 5 (A2) 

where 

&20CO = 0j3(Xiia 5 i + f3 j4 ,a 2 ia 5 i 

- 7/^5 ("21031 + ana4i)/2 , (A3) 
hi(j) = P33 (ana 52 + ai 2 a 5 i) 

+ Pj4 («2ia52 + a22«5l) 

- 7/3j5 [a21«32 + «22Ck31 

+ aiia42 + ai2Q!4i] /2 . (A4) 

Inserting (|38|l . with coefficient a 2 o(j) and an(j) given by 
(|A2I) , into the diagonalized equation (|3"7|) for j — 1 yields 
an expression for the last required coefficient 



&21OO = /Sis [o2o(3) (a 12 a 53 + ai 3 a 52 ) + a 2 o(4) (ai 2 a 54 + a 14 a 52 ) + a 20 (5) (ai 2 a 55 + a 15 a 52 ) + 
an (3) (aiiQ!53 + ai3Q!5i) + an(4) (anas4 + awasi) + an (5) {ana 55 + ai 5 a 5 i)] 
+ ^i4[a2o(3) (a22"53 + a 23 a 52 ) + a 2 o(4) (a22«54 + "240:52) + 
«2o(5) (a 22 a 55 + a 25 a 52 ) + an (3) (a 2 ia 5 3 + a 23 a 5 i) + 
an (4) (a 2 ia 5i + a 2 4a 5 i) + aii(5) (a 2 ia 55 + a 25 a 5 i)] 
- ^/3i5 [O2o(3) (a 22 a 33 + a 23 a 32 + ai2a43 + "13^42) + (A5) 
a 2 o(4) (a 22 a 3 4 + a 2 ^a 32 4- a\ 2 a.^ + ai4«4 2 ) + 
a 20 (5) (a 22 a 35 + ai 2 a i5 + a 25 a 32 + ai 5 a i2 ) + 
an (3) (a 2 ia 33 + aiia 43 + a 23 a 3 i + ai 3 a 4 i) + 
an (4) («2ia34 + 011^44 + 024^31 + 0114041) + 
an(5) (a 2 ia 35 + an 045 + a 25 a 31 + auan)] . 
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